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Abstract. Effect of biquadratic exchange on phase transitions of a planar classical Heisen- 
berg (or XY) ferromagnet on a stacked triangular lattice is investigated by Standard Monte 
Carlo and Histogram Monte Carlo simulations in the region of a bilinear to biquadratic 
exchange interaction ratio J1/J2 ^ 1- The biquadratic exchange is found to cause separate 
second-order phase transitions in a strong biquadratic exchange limit, followed by simulta- 
neous dipole and quadrupole ordering, which is of first order for an intermediate range of 
the exchange ratio and changes to a second-order one again as J1/J2 is further increased. 
Thus, a phase diagram featuring both triple and tricritical points is obtained. Further- 
more, a finite-size scaling analysis is used to calculate the critical indices for both dipole and 
quadrupole kinds of ordering. 
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1 .Introduction 



It has long been recognized that there are materials in which also biquadratic (or generally 
multipolar) interactions may play an important role as in some cases they are comparable 
or even much larger than the bilinear ones. The problem of higher-order interactions in 
systems with Heisenberg symmetry has been tackled in several mean field approximation 
(MFA) studiesi-^, by high-temperature series expansion (HTSE) calculations^, as well as 
within a framework of some other approximative schemes^i^. Those studies have shown 
that such interactions can induce various interesting properties such as tricritical and triple 
points, quadrupole ordering, separate dipole and quadrupole phase transitions etc. Much 
less attention, however, has been paid to this problem on systems with XY spin symmetry. 
Chen et.aU-'^ calculated transition temperatures and the susceptibility critical indices for a 
planar Heisenberg ferromagnet with biquadratic exchange on cubic lattices by the HTSE 
method. However, they only investigated the region of Ji/ J2 ratio for which MFA predicts a 
second-order transition and hence, they failed to provide a complete phase diagram with all 
possible phase transitions and their nature for this model. For example, since MFA predicts 
first-order transitions for dipole ordering if J1/J2 < 1, only quadrupole phase transitions, 
which are expected to be of second order^, were investigated. Here we note that, in spite of 
all the previous conclusions based on various approaches, the rigorous proof of the existence 
of dipole long-range order (DLRO), corresponding to the ferromagnetic directional arrange- 
ment of spins, and quadrupole long-range order (QLRO), representing an axially ordered 
state in which spins can point in either direction along the axis of ordering, at finite temper- 
ature on the classical bilinear-biquadratic exchange model has only recently been provided 
independently by Tanaka and Idogaki^*^, and Campbell and Chayes— . 

We have recently performed MC simulations on the XY model with the bilinear- 
biquadratic exchange Hamiltonian on a simple cubic lattice in order to establish phase 
boundaries and describe critical behaviour of the systen>i^. Unfortunately, we did not suc- 
ceed in finding a conclusive answer to the problem of the order of the DLRO transition in 
the region of the separate DLRO and QLRO transitions. In the present paper we chose 
the lattice with a hexagonal (stacked triangular) symmetry, increased the lattice size and 
performed a careful finite-size scaling (FSS) analysis in order to obtain more reliable data 
and deliver a complete picture of a long-range ordering and its nature for the considered 
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system via Standard Monte Carlo (SMC) and Histogram Monte Carlo (HMC) simulations. 
The phase diagram obtained captures all important features induced by the biquadratic 
exchange such as separate ordering of dipoles and quadrupolcs, the appearance of triple 
and tricritical points, short-range order above the quadrupole ordering temperature, etc. 
Furthermore, we perform a FSS analysis to determine the order of a transition from the 
energy cumulant scahng, and to calculate the critical indices for both DLRO and QLRO 
transitions. 

2. Model and simulation technique 

We are concerned with the planar classical Heisenberg model with bilinear and biquadratic 
interactions, described by the Hamiltonian 

H^-j^j2Si- - J, J2{s, ■ s,r , (1) 

where Si — {Six, Siy) is a two-dimensional unit vector at the ith lattice site , denotes 
the sum over nearest neighbors, and Ji, J2 > are the bilinear and biquadratic exchange 
interaction constants, respectively. 

We first perform SMC simulations on systems oi N = spins, where L = 12, 15, 18, 24, 
30, assuming periodic boundary condition throughout. For a fixed value of the exchange ratio 
Ji/ J2 we run simulations starting at low (high) temperatures from a ferromagnetic (random) 
initial configuration and gradually raise (lower) temperature. These heating-cooling loops 
serve to check possible hysteresis, accompanying first-order transitions. As we move in 
(J1/J2, ksT/ J2) space, we use the last spin configuration as an input for calculation at the 
next point. We sweep through the spins in sequence and updating follows a Metropolis 
dynamics. In the updating process, the new direction of spin in the spin flip is selected 
completely at random, without any limitations by a maximum angle of spin rotation or 
allowed discrete set of resulting angle values. Thermal averages are calculated using at most 
1 X 10^ Monte Carlo steps per spin (MCS/s) after equilibrating over another 0.5 x 10^ MCS/s. 
We calculate the system internal energy E and some other physical quantities defined as 
follows: 

the specific heat per site c 
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the DLRO parameter (magnetization) m, 

m = (M) /N, where M = 
the QLRO parameter q, 

q = {Q)/N, where Q = 



'S'ja, j +1 ^ S'jj^ 



^ (('S'ix) — (Siy) ) 1 + I 2SixSiy 



and the following quantities which are functions of the parameter O 
the susceptibility per site xo 



M, Q) 



Xo 



NkeT 



the logarithmic derivatives of (O) and (O^) with respect to K = l/fc^T 

Dio 



D 



dK 
d 



(O) 
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the fourth-order long-range order (LRO) cumulants Ui and U2 



^1 1 3^02)2 ' 

and the fourth-order energy cumulant V 

V 



Uo 



(3) 



(4) 



(5) 

(6) 
(7) 

(8) 



(9) 



3(^2)2 

These quantities are used to estimate the nature of a transition as well as its position. First- 
order transitions usually manifest themselves by discontinuities in the order parameter and 
energy, and hysteresis when cooling and heating. If transition is second order, it can be 
approximately localized by the xo peak position and also double-checked by the position of 
the fourth-order LRO and energy cumulants curves intersection for various lattice sizes. 

In the next stage we perform HMC calculations, developed by Ferrenberg and 
Swendsen^i^, at the transition temperatures estimated from the SMC calculations for each 
lattice size. Here, 2 x 10^ MCS are used for calculating averages after discarding another 
1 X 10^ MCS for thermalization. We calculate the energy histogram P{E), the order pa- 
rameters histograms P{0) (O = M,Q), as well as the physical quantities (|2])-(l9]). Using 



data from the histograms one can calculate physical quantities at neighboring temperatures, 
and thus determine the values of extrema of various quantities and their locations with high 
precision for each lattice size. In such a way we can obtain quality data for FSS analysis 
which determines the order of the transition and, in the case of a second-order transition, 
it also allows us to extract critical indices. For example, the energy cumulant V exhibits a 
minimum near critical temperature Tc, which achieves the value V* = ^ in the limit L — t- oo 
for a second-order transition, while V* < ^ is expected for a first-order transition^^ii^. The 
extrema of a variety of thermodynamic quantities at a second-order transition are known to 
scale with a lattice size as, for example: 

XO,ma.{L) OC L-^o/uo ^ (10) 

Dio,maAL) L'/"'' , (11) 
^2O,max(^)0cL^/^O , (12) 

where uq and 70 represent the correlation length and susceptibility critical exponents, re- 
spectively. In the case of a first-order transition (except for the order parameters), they 
display a volume-dependent scaling, oc L^. The simulations were performed on the vector 
supercomputer FUJITSU VPP700/56. 

3. FSS analysis and phase diagram 

We started our calculations in a region where the biquadratic exchange is equal to the 
bilinear one, i.e. J1/J2 = 1, and gradually lowered the ratio down to 0. For the ratio greater 
than J1/J2 — 0.55 we found only one phase transition corresponding to a dipole long-range 
ordering which is of second order. Hence, as far as phase transitions are concerned, in this 
range of the exchange ratio the presence of the biquadratic exchange is only reflected in 
the position of the transition temperature (as well as critical indices to some extent) but 
causes no quahtative changes compared with the case of J2 = O^^"^. The second-order 
character of a transition in this region was also confirmed by HMC calculations. Neither 
double-peaked energy or order parameter histogram nor volume-dependent scaling of the 
calculated quantities were observed. Fig.l shows the FSS calculation of z/j\/ from Eqs. ( ITTj) 
and (IT^ . and 7m from Eqn. (ITU]) in a In-ln plot, for J1/J2 = 0.8. The best fit value 
z/M = 0.677 ± 0.008^ agrees very well with the values obtained from other calculations for 
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J2 = 0, but 7Af = 1.177 ± 0.015 is apparently lower than the values obtained for J2 = 
{um = 0.669 and •jm = 1.316 from Refs.— ^i^). However, the critical indices are expected to 
change by the presence of the biquadratic exchang and, as shown by the series-expansion 
calculations for J1/J2 > 1^, 7a/ is continuously lowered by increasing biquadratic exchange. 

As the exchange ratio is lowered slightly below J1/J2 — 0.55, the order of the transition 
changes. Still no separate quadrupole ordering is observed, however, the transition becomes 
apparently first order. This is clearly seen in Fig. 2, where bimodal nature of energy distri- 
bution histograms for various lattice sizes at J1/J2 = 0.5 is obvious. With increasing lattice 
size the barrier between the two energy states is increasing, indicating that the energy is 
discontinuous and hence, the transition of first order. However, as can be seen from Fig. 3, 
for the case of J1/J2 = 0.35 the bimodal distribution can only be observed at sufficiently 
large L, indicating vicinity of the multicritical point. We note that similar behaviour was 
also observed in both order parameters distribution histograms, although we do not show it 
here. 

The bimodal distribution vanishes as the ratio drops below the value Ji/ J2 — 0.33. More- 
over, below J1/J2 — 0.33 quadrupoles start ordering separately at temperatures higher than 
those for dipole ordering. Thus the phase boundary branches and a new middle phase of ax- 
ial quadrupole long-range order (QLRO) without magnetic dipole ordering opens between 
the paramagnetic and DLRO phases. This phase broadens as J1/J2 decreases, since the 
QLRO branch is little sensitive to the J1/J2 ratio variation and levels off, while the DLRO 
branch turns down approaching the point {Ji/ J2-,kBT / J2) = (0,0). This means that the 
ground state is always magnetic as long as there is a finite dipole exchange interaction. In 
Fig.4 we present the temperature variation of both DLRO and QLRO parameters m and 
g, respectively, at J1/J2 = 0.1. We can see that quadrupoles order before dipoles, forming 
a fairly broad region of QLRO without DLRO. The question is of what order are these 
transitions. MFA predicts a first-order transition to DLRO phase in this region^"^, although 
some other MFA results allow possibility of a second-order DLRO transition in the high 
biquadratic exchange limiti. Unfortunately, no other calculations by more reliable tech- 
niques are available. In fact, determination of the order of a transition here turned out to be 
quite tricky. Although no discontinuities in the internal energy temperature variation were 
observed, extremely sharp magnetization slopes and consequently, extremely sharp suscepti- 
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bility peaks at the transition (Fig. 5) would rather suggest a first-order transition. Moreover, 
the fourth-order cumulant U2 at the transition falls to negative values displaying a minimum, 
i.e. displays behaviour seemingly typical for a first-order transition. This, however, must be 
taken with precaution since, here, with increasing temperature the system does not enter 
a paramagnetic phase but a different ordered phase. Indeed, further investigation showed 
that the minimum did not scale with the lattice size as it should in the case of a first-order 
transitions^ and, hence, that in this case it is of different origin and can not be seen as an 
indicator of a first-order transition. Therefore, in order to resolve the ambiguities in the 
phenomena observed from SMC calculations we performed a finite-size scaling analysis from 
HMC data. As mentioned in Introduction, a first-order transition can be reliably detected 
via the fourth-order energy cumulant V scaling, which is shown in Fig. 6 for J1/J2 = 0.1. As 
we can see, a minimum near extrapolated to the limit L ^ 00, V*, reaches the value very 
close to |, as it should be in the case of a second-order transition. Also the other quantities 
clearly did not scale with volume. Instead, as shown in Fig. 7, the obtained scaling indices 
um = 0.642 ± 0.006 and 7m = 1-241 ± 0.011 match quite well the values expected at a 
second-order transition for the given quantities and universality class (See comments on the 
universality class in question in the Discussion and concluding remarks section). Although 
the transition at J1/J2 = 0.1 appears to be undoubtedly of second order, it might change to 
a first-order one at higher ratio J1/J2 but before the DLRO and QLRO boundaries merge, 
as suggested in Ref.-. Therefore, we checked for such a possibility but did not find it to 
be so. Scaling analysis, similar to that for the case of J1/J2 = 0.1, performed for the ratio 
-^^17-^2 = 0.3, which is close to the limiting value J1/J2 — 0.33, produced V* = 0.6666(4), a 
value which deviates slightly more from | than in the previous case, but is still too close to 
be considered as a sign of a first-order transition. Also the critical indices um = 0.654±0.009 
and 7m = 1.216±0.006 confirm the second-order nature of the transition. Therefore, we con- 
clude that the DLRO transition remains second order in the whole region up to J1/J2 — 0.33. 

On the other hand, the nature of the separate QLRO transition is much more obvious. 
All calculated observables behave typically like those for a second-order transition and so 
does the scaling. This is in agreement with the work of Carmesin^, who performed mapping 
of quadrupolar interactions to dipolar ones for planar systems and, hence, showed that a 
transition in these systems should be second order. Critical indices describing a QLRO 
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transition at J1/J2 = 0.2 are extracted from the scaling depicted in Fig. 8. Their magnitudes 
vq = 0.660 ±0.008 and 7q = 1.305 ±0.011 are in a good agreement with the values obtained 
for critical indices connected with a DLRO transition. These values are also in a good agree- 
ment with the values obtained from the high-temperature series expansion calculations for 
cubic lattices^ii^. 

Observing the specific heat temperature dependence we noticed a shoulder-like broad- 
ening just above the quadrupole-ordering transition temperature (Fig.9). We ascribe this 
phenomenon to short-range order, which was already predicted to exist by Chen and Levy^, 
however, apparently could not be seen within their MFA scheme. Here, we show a snapshot 
of the spin arrangement in a XY plane taken at Ti > Tq (inset). In the snapshot we can see 
that the quadrupoles are not totally disordered but tend to align locally in regions equally 
partitioned among the crystalographically equivalent three directions. At the quadrupole- 
ordering temperature the quadrupoles choose one ordering direction and the entire crystal 
forms a single domain. 

The resulting phase diagram is drawn in Fig. 10 and critical indices calculated for second- 
order transitions are summarized in Table [H 

5. Discussion and conclusions 

In this work we studied effects of the biquadratic exchange on the phase diagram of the 
classical XY ferromagnet on a stacked triangular lattice (STL). We believe that our inves- 
tigations, which are the first of this kind for a system with hexagonal lattice symmetry, 
covered all significant phenomena induced by the presence of the biquadratic exchange and 
present a fairly compact picture of the role of this higher-order exchange interaction on 
the critical behaviour of the system considered. We obtained the phase diagram with two 
ordered phases, featuring some interesting phenomena such as the appearance of tricritical 
and triple points. We found that in the region where the bilinear exchange is dominant 
there is only one phase transition to the DLRO phase, which remains second order until 
the exchange ratio reaches the value J1/J2 — 0.55. Upon further lowering of the ratio, 
the transition changes to a first-order one at the tricritical point and remains so down to 
J1/J2 — 0.33. Below this value the phase boundary splits into the QLRO transition line at 
higher temperatures and the DLRO transition line at lower temperatures, which are both 
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of second order. Hence, there is a triple point at J1/J2 — 0.33, where two second-order and 
one first-order transition boundaries meet. 

We speculate that the mechanism driving the system through the mentioned transitions 
could be as follows. At J1/J2 — 00 (J2 = 0) the transition is known to be second order. If 
we introduce the biquadratic exchange and reduce J1/J2, i-e. reduce the bilinear couplings, 
this might induce in a sense a kind of tension (competition) between the two exchange inter- 
actions. Namely, while the decreasing bilinear exchange drives the transition temperature 
down to the lower values, the biquadratic exchange does not follow this tendency and rather 
prevents the ordering temperature from rapid decrease. This tendency is clearly seen from 
the phase diagram both in the region of separate transitions, where Tq does not vary much 
with decreasing J1/J2, as well as in the region of simultaneous ordering, where the transition 
temperature is apparently enhanced by the presence of the biquadratic exchange (the case 
of absent biquadratic exchange is represented by the dash-dot straight line in (Ji — ksTc) 
parameter space). Put differently, quadrupoles would prefer ordering at higher temperatures 
but as long as there is a single transition they are prevented to do so by too low bilinear 
exchange, and order occurs only if the temperature is lowered still further. This frustration 
is enhanced as J1/J2 decreases, and below J1/J2 — 0.55 it results in a first-order transition 
when the strength of the quadrupole ordering prevails and frustrated quadrupoles order 
abruptly along with dipoles. However, below Ji/ J2 — 0.33 the frustration is too high for the 
two kinds of ordering to occur simultaneously - they separate, the tension is released and 
the transitions become second order again. 

The calculated critical indices in the region < J1/J2 < 0.33 for QLRO and 0.55 < Ji/ J2 
for DLRO transitions were found to belong to the standard three-dimensional XY univer- 
sality class, although they showed a slight variation with changing J1/J2, as had also been 
observed in the HTSE calculationa^i^. The values of the indices for the DLRO transition in 
the region < J1/J2 < 0.33 appear to be in a reasonable agreement with those just men- 
tioned above, and one might tend to include them into the same XY universality class. We 
believe, however, that they rather belong to the Ising universality class (Note that the critical 
indices of both universality classes take fairly similar values: u^"^ = 0.669, 7"^^ = 1.316 and 
j^ismg _ Q_g29, 'j^'^'^^d = 1.239^^ for J2 = 0). The reason for our claim is that in this region 
there is first axial quadrupole ordering and only then directional dipole ordering within the 
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given axis. Therefore, the only difference from the Ising case is that dipoles can order along 
any of the three axes, not only the z-axis. The Ising-like nature of the ordering would also 
explain the remarkably sharp behaviour of the physical quantities at the transition as, for 
example, shown in Fig.5. 

We further intend to perform similar simulations on STL with antiferromagnetic bilin- 
ear and antiferromagnetic/ferromagnetic biquadratic exchange. Such a spin system would 
present a geometrically frustrated system with additional frustration arising from compe- 
tition between the antiferromagnetic bilinear and ferromagnetic biquadratic, or inter-plane 
antiferromagnetic bilinear and antiferromagnetic biquadratic exchange interactions. Among 
the main points of interest of such investigations will be the nature of the ground-state order, 
the universality class of the antiquadrupolar ordering, etc. 
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FIG. 1: Scaling behaviour of the maxima of the susceptibihty XM,max and logarithmic derivatives 
of the DLRO parameter and its second moment Dim, max and D2M,max, respectively, in In-ln plot, 
for J1/J2 = 0.8. The slopes yield values of 1/um for Dim, max, D2M,max and ^m/vm for XM,max- 

TABLE I: Critical indices vq^'^q] vnuIm] and z/, 7 for quadrupole, dipole and simultaneous order- 
ing, respectively. 
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0.654 ± 0.009 
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FIG. 2: Energy distribution at the size-dependent transition temperatures Tc{L) for various lattice 
sizes and J1/J2 = 0.5. Double-peaked structure with deepening barrier between the two energy 
states with increasing lattice size indicates a first-order transition. 
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FIG. 3: Energy distribution at the size-dependent transition temperatures Tc{L) for various lattice 
sizes and J1/J2 = 0.35. The bimodal distribution signahng a first-order transition can only be seen 
at sufficiently large L. 
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FIG. 4: Temperature variation of the DLRO and QLRO parameters m and q, respectively, for 
J1/J2 = 0.1 and L = 12. 
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FIG. 5: Susceptibility xm vs temperature for different values of J1/J2 and L = 18. 
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FIG. 6: Scaling of the energy cumulant minima with volume at J1/J2 = 0.1. The value V* = 
0.66666(5) is obtained from the extrapolation L — ^ 00. We note that the brackets here only signify 
that the last digit could not be determined with certainty from the best fit and, hence, the value 
in the brackets should not be seen as an error estimate. 
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FIG. 8: Scaling behaviour of the maxima of the susceptibihty XQ,max 

and logarithmic derivatives 

of the QLRO parameter and its second moment Diq-max and D2Q^max, respectively, in In-ln plot, 
for J1/J2 = 0.2. The slopes yield values of I/vq for Dig^rnax, ^2Q,max and 'Jq/vq for XQ,max- 
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FIG. 9: (a) Specific heat vs temperature showing a shoulder-like broadening at T\ >Tq, indicating 
the presence of short-range order above a quadrupole ordering temperature, (b) Snapshot shows 
signs of a local alignment of spins along the crystalographically equivalent three directions. 
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FIG. 10: Phase diagram in (J1/J2) fcs7c/'^2) space. The solid and dashed Hnes correspond to 
second- and first-order transitions, respectively, and the dash-dot straight hne represents the bound- 
ary between paramagnetic and ordered regions in (Ji, /cgTc) space when the biquadratic exchange 
is absent. 
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